library(dplyr)
library(lme4)
library(moments)
library(modelsummary)
library(broom.mixed)

data<-readRDS("table4_data.rds")


services <- c("Waste_treatment", "Waste_collection", "Transport","Fire",  "Libraries",
  "Civil_protection",  "Drinking_water", "Sewer")
 

transform_skewed_data <- function(data) {
  # Function to handle skewness transformation
  transform_variable <- function(x) {
    skew <- skewness(x, na.rm = TRUE)
    
    if (abs(skew) > 1) {
      if (skew > 0) {
        # For positive skew, try log transformation if all values are positive
        if (min(x, na.rm = TRUE) > 0) {
          x_transformed <- log(x)
        } else {
          # If there are zero or negative values, shift the data
          x_shifted <- x - min(x, na.rm = TRUE) + 1
          x_transformed <- log(x_shifted)
        }
      } else {
        # For negative skew, try exponential transformation
        x_transformed <- exp(scale(x))
      }
      return(as.vector(scale(x_transformed)))
    } else {
      return(as.vector(scale(x)))
    }
  }
  
  # Apply transformations and scaling
  data %>%
    mutate(
      population = transform_variable(population),
      population_squared = transform_variable(I(population^2)),
      VOTANTS_PERCENT = transform_variable(VOTANTS_PERCENT),
      EUROS_HABITANT = transform_variable(EUROS_HABITANT),
      year = factor(year)
    )
}

# Modified model fitting function to include both random effects
safe_glmer_both <- function(service_name, data) {
  service_data <- data %>% filter(service == service_name)
  
  tryCatch({
    model <- glmer(imc_dummy ~ population + I(population^2) + VOTANTS_PERCENT + EUROS_HABITANT + 
                     (1 | year) + (1 | muni),
                   data = service_data,
                   family = binomial,
                   control = glmerControl(
                     optimizer = "bobyqa",
                     optCtrl = list(maxfun = 2e4),
                     check.conv.grad = .makeCC("warning", 0.2),
                     check.conv.singular = .makeCC(action = "warning", tol = 1e-4)
                   ))
    
    if(isSingular(model)) {
      cat("Singular fit for", service_name, "\n")
    }
    
    return(model)
  }, error = function(e) {
    cat("Error in", service_name, "model:\n")
    print(e)
    return(NULL)
  })
}


# Transform the data with skewness handling
data_transformed <- transform_skewed_data(data)

# Alternative approach using only base R
service_names <- unique(data_transformed$service)
models_list_both <- setNames(
  lapply(service_names, function(x) safe_glmer_both(x, data_transformed)),
  service_names
)


model_summaries <- models_list_both[!sapply(models_list_both, is.null)]

# Function to check singularity and convergence issues
check_model_issues <- function(models) {
  issues <- lapply(names(models), function(name) {
    model <- models[[name]]
    if (!is.null(model)) {
      singular <- isSingular(model)
      conv <- model@optinfo$conv$lme4$messages
      c(Singular = singular,
        Convergence = !is.null(conv),
        Messages = if(!is.null(conv)) paste(conv, collapse="; ") else "None")
    }
  })
  names(issues) <- names(models)
  return(do.call(rbind, issues))
}

# Get model diagnostics
model_issues <- check_model_issues(model_summaries)
print("Model Fitting Issues:")
print(model_issues)

# Create summary table
modelsummary(
  model_summaries,
  stars = TRUE)

# Save results to Word document
modelsummary(
  model_summaries,
  stars = TRUE,
  output = "mixed_effects_submitted.docx"
)